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Abstract 

We introduce a general Hamiltonian describing coherent superpositions of Cooper pairs and 
condensed molecular bosons. For particular choices of the coupling parameters, the model 
is integrable. One integrable manifold, as well as the Bethe ansatz solution, was found by 
Dukelsky et al, Phys. Rev. Lett. 93 (2004) 050403. Here we show that there is a second 
integrable manifold, established using the boundary Quantum Inverse Scattering Method. In 
this manner we obtain the exact solution by means of the algebraic Bethe ansatz. In the 
case where the Cooper pair energies are degenerate we examine the relationship between the 
spectrum of these integrable Hamiltonians and the quasi-exactly solvable spectrum of particular 
Schrodinger operators. For the solution we derive here the potential of the Schrodinger operator 
is given in terms of hyperbolic functions. For the solution derived by Dukelsky et al, loc. cit. the 
potential is sextic and the wavefunctions obey T^T-symmetric boundary conditions. This latter 
case provides a novel example of an integrable Hermitian Hamiltonian acting on a Fock space 
whose states map into a Hilbert space of "PT-symmetric wavefunctions defined on a contour in 
the complex plane. 



1 Introduction 

One of the most interesting developments in the physics of ultra-cold gases has been the production 
of molecular condensates from fermionic constituents [1,2]. Due to the property of Pauli blocking, 
it turns out that molecular condensates formed from fermionic atoms are generally more resistant 
to decay compared to the case of molecules formed from bosonic atoms. Furthermore, such systems 
open possibilities to experimentally probe the nature of the crossover between BCS and BEG type 
physics [3]. 

In order to theoretically model such systems, it is useful to have integrable systems with exact 
solutions which allow for investigations beyond the limits of perturbative and mean-field treatments. 
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Recently, Dukelsky et al. [4, 5] have proposed such a model and derived an exact solution using the 
methods of Richardson-Gaudin type spin models [6]. In their approach, a bosonic degree of freedom 
describing the molecular condensate is introduced by taking the infinite spin limit of one realisation 
of particular spin operators following a procedure of Gaudin [7]. The remaining spin operators are 
then realised in terms of paired creation and annihilation operators, which may be either bosonic or 
fermionic. 

The aim of the present work is to derive a similar, although significantly distinct, model which 
also admits an exact solution. The approach we adopt is to use the boundary Quantum Inverse 
Scattering Method (QISM) as developed by Sklyanin [8]. We construct a doubled monodromy 
matrix from the spin realisation of the Yang-Baxter algebra in the usual way but with a bosonic 
operator valued solution of the reflection equations, which is obtained by dressing a boundary K- 
matrix with the bosonic realisation of the Yang-Baxter algebra. For simphcity, here we only consider 
the case where the spin operators are realised in terms of Cooper pairs of spin 1/2 particles, although 
generalisations to higher spin particles are possible. Through this construction we yield a family of 
commuting transfer matrices which give rise to an abstract intcgrablc system. Following the method 
used to establish the integrability of the Russian Doll BCS model [9], we obtain the conserved 
charges by expanding the transfer matrix in inverse powers of the spectral parameter. Finally, by 
taking the quasi-classical limit (e.g. see [10]) we obtain a Hamiltonian that is expressible in terms of 
the conserved charges, thus establishing integrability. We also undertake the algebraic Bethe ansatz 
to derive the associated exact solution. 

The integrable Hamiltonian that we will derive, and a particular case of the ones given in [4,5], 
both belong to the following class of Hamiltonians: 



where the Sj denote M two-fold degenerate energy levels for a system of fcrmions, uj is the single 
energy level for condensed molecular bosons, U is the interaction energy associated with S-wave 
scattering of the molecular bosons, and Qj are the interaction strengths for scattering between Cooper 
pairs and molecular bosons. In this article, the spin operators represent Cooper pairs through a 



Further, 6, W arc canonical Bose operators obeying [6, 6^] = 1. We adopt the notation TV^ = }}^h, 
rij = Uj^ + nj+, Nj = rij^rij-. It can be verified that the total "bosonic" particle number (i.e. the 

M 

number of Cooper pairs plus the number of molecular bosons) N — Nh + Nj commutes with the 
Hamiltonian and is thus conserved. We note that for any j the change of variable 



M 



M 




(1) 




-9j 



(2) 



is a unitary transformation which is induced by the automorphism 
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Similar to BCS models [6, 9, 10] the Hamiltonian displays Pauli blocking. For any unpaired 
fermion at energy level Sj the action of the interactions with coupling gj is zero. This means that the 
Hilbert space can be decoupled into a tensor product of paired and unpaired fermion states for which 
the action of the Hamiltonian on the subsystem of unpaired fermions is automatically diagonal in 
the natural Fock basis. Below we let 1 denote the index set for the unblocked levels and let £ < M 
denote the cardinality of 1. 
It was shown in [4, 5] that 

f/ = 0, g,=g Vj (3) 
is an integrable manifold for the Hamiltonian (P). For this case the energy eigenvalues are given by 

N 

where the parameters Xj are solutions of the Bethe ansatz equations 



E ~ (^) 



29 2 ' 4 2s,/ g - x, - x, 

In the first part of the paper we will show that there is a second manifold of integrability which holds 
when 

and we will derive the corresponding Bethe ansatz solution. We remark at this point that there are 
two types of elementary excitations for the model, those for which Cooper pairs are broken and those 
for which they are not. Clearly for a non-pair breaking excitation C — 2N remains constant. For 
a pair breaking excitation N decreases by one and C decreases by two, such that C — 2N remains 
constant. Hence the coupling constraints defined by (0) are independent of the type of excitation. 

In the second part of the paper we will show that for particular submanifolds of both (jH)) and 
(0) there exists in each case a mapping of the energy spectrum into part of the spectrum of a quasi- 
exactly solvable (QES) Schrodinger equation. Quasi-exactly solvable problems are characterised by 
having a finite number of wave functions and eigenvalues that can be found algebraically [11,12]. 
The submanifolds of (jH)) and (0) we consider map onto the QES sectors of two Schrodinger equations. 
The case of is particularly interesting because it naturally maps onto a PT-symmetric [13, 14] 
Schrodinger problem. 



2 Boundary QISM and construction of the Hamiltonian 

We will construct a family of commuting transfer matrices t{u) via the boundary QISM [8] from 
which the Hamiltonian will be obtained. For the operator t{u) acting on a Hilbert space of physical 
states, we require 

[t{u),t{v)] = 0, W,veC. 
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A consequence of constructing the transfer matrix in this way is that by taking the series expansion 

oo 

tiu) = ^ tkl 

we have 



k=—oo 



which represent constants of the motion. Below we demonstrate how the Hamiltonian is expressible 
as a function of these constants of motion and as a consequence deduce that the resulting model is 
integrable. 

We begin with the Yang-Baxter equation 

Ri2{u — v)Riz{u — w)R23{v - = R23{v — 'w)Ri3{u — w)Ri2{u — v) 
and use the su{2) invariant i?-matrix solution which has the form 



R{u) = 



/I 0\ 

b{u) c{u) 

c{u) b{u) 

\0 1/ 



with b{u) — u/u + T], c{u) = ri/u + T] where rj is an arbitrary complex parameter. We require an 
L-operator, a realisation of the Yang-Baxter algebra acting on local physical spaces, which here is 
given by a spin 1/2 realisation of the su{2) Lie algebra 



for 



We note that for 



we have 



^^J 2 



u 



1 - TV- 

2 ^^J 



(6) 



LJu) — I Sj, 

■' u-rj 



Lj{u)Lj{u) 
The L-operator so defined satisfies 



4:U{U — Tj) 



I. 



Ri2{u — v)Li{u)L2{v) = L2{v)Li{u)Ri2{u — v). 



(7) 



(8) 



We also take the following realisation of the Yang-Baxter algebra in terms of the molecular boson 
operators &, 



J{u) = 
J{u) = 



l + r]u + vi'^{Ni, + l) rjb 
r]b^ 1 

1 -r]b 
-rjb^ l + r]u + rfNb 



where we also have 

J{u)J{u) = {1 + 7]U)I. (9) 
The operator J{u) obeys the following relation 

Rl2{u - v)Ji{u)J2{v) = J2iv)Ji{u)Ri2{u - v) . (10) 

For further details on these realisations of the Yang-Baxter algebra, the reader is referred to [10]. It 
will be useful to define 

We construct the doubled monodromy matrix [8] in terms of the local operators for Cooper pairs 
and molecular bosons. For a given index set X for the unblocked levels with cardinality C, we first 
perform a relabelling such that I ^ I' = {1, 2, — 1, £}. The doubled monodromy matrix is 

given by 

T{u) = Li{u-ei)...Lc{u-ec)J{u)KJ{-u)Lc{-tc-u)...U{-ti-u) (12) 

where 



K 



1 
-1 



is the boundary i^'-matrix. Note that here we will choose the boundary i^'-matrices, usually denoted 
Kj^, K^, to be equal. These are not the most general ii'-matrices that can be chosen, but this choice 
is sufficient for the considerations of the present work. The i^-matrix satisfies the equation 

Ru{u -v){K® I)Ri2{u + v){I(^K) = {K® I)Ri2{u + v){K ® I)Ru{u - v). (13) 

The contribution J{u)KJ{—u) to the doubled monodromy matrix is a dressing of the boundary 
ii'-matrix, which itself can be considered as an operator valued boundary ii'-matrix. We remark 
that operator valued boundary if-matrices have been used previously for the purpose of embedding 
Kondo impurities in integrable one-dimensional t — J models [15,16]. 

From (j7l8l9ll0ll3|) it follows that the doubled monodromy matrix satisfies 

Ri2{u - v){T{u) ® I)Ri2{u + v){I® T{v)) = (/ ® T{v))Ri2{u + v){T{u) ® I)Ri2{u - v). (14) 
Taking the expansion of the doubled monodromy matrix in inverse powers of u, we can express it as 

T{u) ^ uTo + u^Ti + u-^T2 + u-^T3 

where 

To = r// 

c 

j 

c c c / c c c 

T2 = r/J](5,ir + ir5',) + 2r72 5^(5,X + X5,)-r73^5, + r7M + + 

i i i \j<i j>i i,j 



c c / c c c \ 

i i \j<i j>i i,j / 



The transfer matrix is given by 

t{u) = Tt[KT{u)] 
which provides a family of commuting matrices: 

[t{u),t{v)] = \/u,v. 
A series expansion of the transfer matrix provides a set of mutually commuting operators, 

t{u) ^ Uto + + U~^t2 + U~^ts, 

[tk,t,]=0. 

We find 
to = 

c 

h = 2 + 2r/2(2iVb + l) + 2r72^(2Ar, -1) 

j 

t2 = 

c c c 

h = VX]e,(5r6t + 5;6)+2r/2X]e,2(2iV,-l)+r/2^ 



(15) 



2{N,--m--) + sjst + sls; 



C- r 



mj-^m-\)+sjs*+s*s: 



WE 

The Hamiltonian is related to by 



C r 



2{N,--m--)-SrSt-S^S, 



i7 = lim ^ (^3 + 2,^ ± el - A.^N^N - C) - 



(16) 



after making the substitutions 



, = Uei, u = U{C- 2N) 



(17) 



where z G X' = {1, . . . £} label the unblocked levels. The above conditions may be easily relabelled 
in terms of the original index set X. 
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3 Algebraic Bethe ansatz solution 

Sklyanin's boundary QISM for the doubled monodromy matrix [8] was adapted for the Gaudin 
magnet in [17, 18] to obtain a solution using the algebraic Bethe ansatz (ABA) method. Using this 
general scheme, we outline the boundary ABA to develop a solution for the Hamiltonian subject 
to the constraints 

With the above definition for the transfer matrix p5|) . we solve the eigenvalue problem 

= A(n)^ 

by the boundary ABA method. Representing the doubled monodromy matrix as 



Tin) 

the transfer matrix is given by 



A{u) B{u) 
C{u) D{u) 



t{u)^ = {A{u)- D{u))^. 

The relations resulting from (fT^ are quite complicated in terms of which contributions are relevant 
to the diagonalisation of the transfer matrix. For this reason, we define a new A{u) through 

A{u) = {2u + rj)A{u) - rjD{u). 

We can now rewrite the relevant relations arising from (|14j) in a more convenient form: 

A( \n( \ {u - V + rii){u + V + 2rii) ^ 2r/(M + r/) ^ 

+ , '"^,1^ + ^^ M u)D{v) (19) 
(2t> + ri)[u + V + r]) 

D{u)C{v) = (; + "^(\;-^| c(.)D(^) + , ^-^ -Mu)D{v) 

^ ' ^ ' (^u + v + ri){u-v) ^ ' ^ ' l^2v + ri){u - v) ^ ' ^ ' 

C{u)A{v). (20) 



{u + V + ri){2v + 7]) 

Correspondingly, the transfer matrix may now be expressed as 

1 ^ 2v 

t(u) = A(u) —D(u). 

^ ^ 2u + r] ^ ^ 2u + ri ^ ^ 



For the construction of the eigenstates, the reference state is chosen to be 

|fi) = |0)(g)(j) (21) 



i=l 



such that B{u)\Q) = 0. We find that the action of A{u), D{u) on this reference state is 

A{u)\n) = 2{u + ri)a{u)a{-u)\n) (22) 
D{u)\Q) = -6{u)6{-u)\Q) (23) 
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where 



i=l ^ 



i=l 

C 



i=i ^ 



6{-u)\n) = {l-vu)\ \ [l + wr^]\^)- 



2(u + ti 



Taking the Bethe ansatz states to be a product of the creation operators on the reference state, 

AT 



we find that for 

= A(n)\I/ + "unwanted terms" 

the co-efficient A(m) is given by 

^^^^ ^ 2{u + r]){l + r)u + T]'^) -pj {u - + r]){u + Vg + 2r]) {u-Ej- r]/2){u + + 7]/2) 
2u + r] {u + Va + ri)iu-Va) {u- ei){u + ei + T]) 

Cancellation of the unwanted terms above, which is equivalent to imposing that A{u) has no poles, 
leads to the Bethe ansatz equations 

C N 

vp{l - r]Vf}) Y\ + g» + 3?7/2)(t;/3 - + 7]/2) ^ -|-|- {vp -Vq + f]) {vp + Vg + 2r]) 
{vp + r]){l+7]V(} + r]^)fJ^{vf3 + ei + r]/2){vp~ei-r]/2) ^/p^^i (vfs + Va)ivp - Va - v) 

(3 = 1...N. (25) 

3.1 Results in the quasi- classical limit 

Expanding the eigenvalue of the transfer matrix (j24j) in inverse powers of u we obtain 



A(u) ^ -2r]^jZ + Arj^N + 2 + 2/7^ + ^ 



£2 - 4£iV - 2 ^ + 47V2 _ + 4 ^ vl 



i=l a=l 
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and in particular 



N 



3C 



a=l 



i=l 



In view of this leads us to the energy of the Hamiltonian being 



N 



k(f_J 



where X again denotes the index set for the unblocked levels. The corresponding Bethe ansatz 
equations arising from (j25p in the quasi-classical limit are 



1 >^ 1 



N 



P fcex 



Y f 



(26) 



4 Mapping to a Schrodinger equation 

For the remainder of the paper we will restrict to a degenerate case for which Sj is independent of 
j and consider only non-pair breaking excitations such that C = M. For both exact solutions given 
by the Bethe ansatz equations (BAE) ()4l26p . we will demonstrate that the spectrum of the model, 
with a suitable scaling, can be mapped into that of a one-dimensional Schrodinger equation. 



4.1 A first mapping 



We begin with the BAE (I26j) derived in the previous section. First observe through (jl7|) that the 
variable U determines the overall scaling of the model. We set U = — 1 (so the S-wave scattering is 
attractive) and also 



The BAE are now 



2ei = 7,Vj, 



1 M 
1 - — + 



2xis xp-'y 



N 



It is convenient to shift the roots, xp ^ Xj3 + 7/2. Then the BAE become 

N 



- 1 



1 



+ 



M 



2{x^ + ^/2) x^-7/2 



^Xp-Xc 



(27) 



Following the procedures taken in [19,20], we choose the following ansatz for an ordinary differential 
equation (ODE): 



F{y) = [y' - {il2f]Q"{y) 



--{y - 7/2) + M{y + 7/2) - y2 + (^/2)^ 



Q\y)- 



{2i 
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N 

For Q{y) = Y\_(y ~ ^p)-> ^^"^ check that F{xp) = using the BAE (P7j) and noting 

/3=1 



^ 1 

Q'iy) = Qiy)Y.7^ 



N 

y -^a 

1 



Q = l 

N N 



The degree of the polynomial Q{y) is N, so in the RHS of (j^H|l the highest order term is + 1. Thus 
we can fix F{y) to be 

F{y) = {Ay + B)Qiy). 
It can be shown that the leading order term gives 

A = N. 

Next consider 

Q(7/2) ^ - 7/2 

Q(-7/2) ^^« + 7/2' 

From the ODE (|^. we also see that 

Q^(7/2) ^ 

Q(7/2) 2M M7 

QX-7/2) ^ ^_2B 
Q(-7/2) 7 ■ 



Rearranging for B we find 

5 



M7Q'(7/2) 7Q'(-7/2) 



2 g(7/2) 4g(-7/2) 

iV27 X ^ 1 7 X ^ 1 

"^^^/3-7/2'^ 4^x^ + 7/2' 

We can simplify B using an identity from the BAE. Taking the BAE (jTTj) . multiplying by a;^ and 
taking the sum over /? gives us 

N N 

-g., + (-l/2.M)A^.2g^.^g^^;V(A^-l). ,29) 
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Using this identity, we obtain B as follows 



B 



-E ~ N M - N ■ 



iV7 



where for the energy we have 



N 



N 



13=1 



/3=1 



Hence Q{y) satisfies the differential equation 
[Ny-E - 



--{y- 7/2) + M{y + 7/2) - + (7/2) 



Q{y) = iy' - {i/'2f)Q"{y) 
Q'{y). 



(30) 



Now we look to put the ODE (j3(jp into the form of a Schrodinger equation. We start with the 
substitution y = — (7coshx)/2, so that 



dQ ^ ■ x^f \dQ 



Now ()30|) becomes 



dy'^ 72sinh^(x) ydx"^ 
d'^Q sinh(x) 



— coth(xj— — 
dx 



dQ 



J 9 0/ t./ A ,T[l + 2M + 7(cosh(x) + l)] ^ 
dx^ 2(cosh(x) + 1) ^ dx 



'^(cosh(a;) \\)-E-N\M-N\- 



Q(x). 



Next we put 

\E'(x) = exp(— (7(x))(5(a;) 

and we need to find g{x) such that upon substituting \I'(x) into the ODE (j28|l . the \E''(x) term cancels. 
(Note that the '^[x) here is not related to the state vector \l/ introduced earlier in the ABA section.) 
For 

Q!{x) = exp{g{x))[g'{x)'^{x) + '^'{x)] 

Q"{x) = exp{gix))[{g"ix)+g'{xf)^{x) + 2g'ix)^'ix) + ^"ix)], 
we find that the contribution to \E''(a;) (which we need to eliminate) is as follows: 

sinh(x) 



This is easily solved for g{x): 



4(cosh(a;) + 1) 



[1 + 2M + 7(cosh(x) + 1)] . 



g[x) = -[7cosh(a;) + (2M + 1) ln(cosh(x) + 1)]. 
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We now obtain 



^(x) = (cosh (x) + i)-(A^/2+i/4) exp (^-1 cosh(x)) Q{y) 



and also note 

"( _ 7cosh(x)(l + cosh(x)) + 2M + 1 
^ ^"^^ ~ 4(cosh(x) + 1) ■ 

Next we can rewrite the ODE (j^H)) as a Schrodinger equation 



where 



^"(x) + V(x)^(x) = E^(x) 



2' ' ' ' 2 J 16 8 

(2M + 1)(2M + 3) ^ (2M - 1)7 cosh(x) ^ 7^ sinh^(x) 



and 



8(cosh(x) + 1) 8 16 



N 

^(x) = (cosh (x) + i)-W2+i/4) exp (^-1 cosh(x)) JJ Q coshx + x„) 



This potential belongs to a class of known quasi-exactly solvable potentials [11,21]. 
4.2 A second mapping 

We can perform a similar mapping for the case of the solution (0}. Setting Sj = e, g 
performing a shift in the roots Xj Xj + e, the BAE are 

M 



2e — u — Xj -\ = — 



where the energy is given by 



N 

E = -^Xj - 2eN. 

We choose the following ansatz for the ODE: 

Fiy) = yQ"{y) + [{co - 2e)y + y^-M] Q\y). 

N 

For Q{y) = Y\(y ~ ^i)' easily check that F{xj) = using the BAE. Now set 

Fiy) = {Ay + B)Q{y). 
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It can be shown that the leading order term gives 

A = N. 



Next consider 



Q(o) 

Q"{^) 
Q(o) 



TV 



N N 



' ^ ' ^ 1^' ■ J 

j=l k<j,k=l ■' 



Also 



^ 1 

F(0) = B = -MQ'{0) = mJ2—- 



We simplify B using an identity from the BAE. Taking the summation over j in the BAE, we obtain 
B as follows 

B = Nlu-E. 
Hence Q{y) satisfies the differential equation 

yQ"{y) + [{u - 2e)y + y' - M] Q\y) + [-Ny ~ Nu + E] Q{y) = 0. (36) 

We make the substitution y = so that 



dQ X dQ 
dx 2 dy 
d^Q 



dy"^ 



A dQ 4 d^Q 
dx dx'^ 



x-^ 



Now the ODE ()36|) becomes 

Nx'^ 



-Nuj + E 



x^ 2M+1 {uj-2e)x 



X 



dQ 
dx 



Next we put 



= exp{—g{x))Q{x) 

and substituting into the ODE, the contribution to \E''(x) is given by 



2g\x) 

which we solve for g{x) to obtain 



x^ (2M +1) (cj - 2e)x 



+ 



X 



, , (2M + 1), x^ x^{uj-2t) 

q(x) = mx 7 . 

^ 2 64 8 
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Setting 2; = x/2, /5 = — 2e, we can now rewrite the ODE ()36|) as a Schrodinger equation 



where 



^"{z) + V{z)^{z) = E^{z), (37) 



Viz) = z' + 2(3z^ + {(3' + 2(1 + 2N-M)}z'+ ^^^ + ^^^^ (38) 

N 



E = -2I3{2N - M) -A'Y^Xj, 



= ^-(^^+^)/^exp('^ + ^m(^^-^a) (39) 



a=l 



and the parameters roots of the BAE 

The quasi-exact solvabihty of the potential (jnH|) has previously been discussed in [12] (see Chapter 
2.2). However some care is needed in order to properly embed the eigenfunctions into a Hilbert space 
of states. Specifically, the difference between the standard Schrodinger problem and the one above 
lies in the boundary conditions imposed on the eigenfunctions in order to obtain the Hilbert space. 
In the standard case the wave functions are required to be square-integrable, and thus decay at ±00 
along the real axis. The wave function (jSHl) is not square-integrable on the real axis, since at large real 
1^1 it clearly blows up. However it is square-integrable if we instead define it on a contour which lies in 
wedges of the complex plane of open angle 7r/4 centered about aigz = —ir/A and arg^; = — 37r/4, and 
distorted to pass below the origin since (M+3/2)(M-|-l/2) 7^ [26]. These are PT-symmetric [13,14] 
boundary conditions, and our QES potential (jHHj) leads to the angular-momentum generalisation of 
the "PT-symmetric Schrodinger equation studied in [27]. It is interesting to note that for (jHHj) . the 
QES spectrum is necessarily real, since it has been derived from a mapping of a particular integrable 
manifold of the Hermitian Hamiltonian (Q). By contrast, the potential studied in [27] admits a region 
of broken PT-symmetry where the QES spectrum becomes complex. We always map to the region 
of unbroken PT-symmetry because of the constraint M > 0. 

Finally, we remark that the QES potential ()38|) has an interesting duality property under the 
transformation (3 —(3, z iz. The quasi-exactly solvable sector of the spectrum maps onto itself 
via E — > —E [23, 24] while the non-QES sectors remain unrelated. In our picture this duality is a 
trivial consequence of the freedom we have to redefine the Bose operators 6 — >■ — 6 and 6^ —bK If 
we send uj —uj and e — > e the Hamiltonian becomes —H and the energy eigenvalues change sign. 
Hence the QES energy levels become —E. 



5 Conclusion 

The general Hamiltonian has two distinct integrable manifolds in parameter space which are given 
by ()3l5p and exact Bethe ansatz solutions given by ()4l26p respectively. For a certain submanifold 
in each instance, we have shown that the eigenspectrum and eigenstates can be mapped to the 
algebraic sector of a QES Schrodinger potential. It is surprising that for one solution the mapping 
is to PT-symmetric eigenstates defined on a contour in the complex plane, while for the other case 



14 



the eigenstates are defined on the real axis. The imphcations of this curious result in relation to the 
difTerent physical properties of the Hamiltonians warrant further investigation. 

Acknowledgements: We thank the Australian Research Council for financial support. 
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